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Abstract 

The present paper points out to a novel scenario for formation of chaotic attractors in a class of models 
of excitable cell membranes near an Andronov-Hopf bifurcation (AHB). The mechanism underlying 
chaotic dynamics admits a simple and visual description in terms of the families of one-dimensional 
first-return maps, which are constructed using the combination of asymptotic and numerical techniques. 
The bifurcation structure of the continuous system (specifically, the proximity to a degenerate AHB) 
endows the Poincare map with distinct qualitative features such as unimodality and the presence of 
the boundary layer, where the map is strongly expanding. This structure of the map in turn explains 
the bifurcation scenarios in the continuous system including chaotic mixed-mode oscillations near the 
border between the regions of sub- and supercritical AHB. The proposed mechanism yields the statistical 
properties of the mixed-mode oscillations in this regime. The statistics predicted by the analysis of the 
Poincare map and those observed in the numerical experiments of the continuous system show a very 
good agreement. 
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Identifying bifurcation scenarios leading to the formation of chaotic attractors is one of the cen- 
tral problems of nonlinear science. In applications, such information can be used for locating and 
characterizing the regions of complex behavior. The present paper identifies the proximity to a de- 
generate Andronov-Hopf bifurcation as a source of complex dynamics in differential equation models 
of Hodgkin-Huxley type. The latter constitute an important class of models of computational biology. 
The results of this work explain the origin of chaotic mixed-mode oscillations generated by models 
of this class and yield the statistical properties of the irregular oscillatory patterns. The statistics 
predicted by the analysis show a very good agreement with those estimated in the numerical experi- 
ments. 
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1 Introduction 



Understanding mechanisms responsible for generating different patterns of electrical activity in neurons, 
firing patterns, and transitions between them is fundamental for understanding how the nervous system pro- 
cesses information. After a classical series of papers by Alan Hodgkin and Andrew Huxley ifTOl , nonlinear 
differential equations became a common framework for modeling electrical activity in neural cells. Today 
the language and techniques of the dynamical systems theory, especially those of the bifurcation theory, 
are an indispensable part of understanding computational biology. A distinctive feature of many biological 
models is that they are often close to a bifurcation. In particular, all excitable systems, including many 
models of neural cells reside near a bifurcation. The type of the bifurcation involved in the mechanism of 
action potential generation underlies a fundamental classification of neuronal excitability ifTTl . According 
to this classification, the excitation in type II models is realized via an AHB. The proximity of a model to 
an AHB can have a significant impact on the firing patterns that they produce. Near the bifurcation, systems 
acquire greater flexibility in generating dynamical patterns varying in form and frequency. The proximity to 
the instability also makes this class of systems more likely to exhibit chaotic behavior. Typical mechanisms 
of transitions to chaotic dynamics in conductance-based models of neural cells are much less understood 
compared to those corresponding to periodic activity. Understanding these mechanisms is important in view 
of the prevalence of irregular dynamics in both modeling and experimental studies of neural cells. 

In the present paper, we describe a regime of chaotic mixed-mode oscillations (MMOs) in type II con- 
ductance based models. Mixed-mode oscillatory patterns combine small amplitude (subthreshold) oscilla- 
tions with large amplitude spikes (Fig. [I]). Despite subtle underlying dynamical structures, MMOs have 
been reported in many experimental and modeling studies of neural cells (3j |U \5\ Q21 [I9j |20l as well as in 
the models of phenomena outside of neuroscience [H [15] EZ1. The mathematical mechanisms generating 
MMOs have recently become an area of intense research. Several structures have been proposed to serve 
as organizing centers for MMOs. Among them, there are folded singularities ID 01231, proximity to the 
AHB combined with the strongly contracting return mechanism lfl"3l . and the Hopf-homoclinic bifurcation 
j9l . Recent studies have greatly elucidated the sources for periodic mixed-mode regimes. The mechanisms 
for irregular MMOs remain to be explained. The present paper suggests a general mechanism for chaotic 
MMOs in terms of the bifurcation structure of the corresponding class of models. Specifically, we show 
that the transition from supercritical to subcritical AHB in a class of systems lies through the regime of 
chaotic MMOs. Consequently, systems lying close to such transition are likely to exhibit chaotic behavior. 
This effect was first noticed in the context of the model of solid fuel combustion in (both in the original 
infinite dimensional model and its finite dimensional approximation) O. The qualitative explanation for 
the mechanism of chaos near the border of criticality (i.e., where the AHB changes its type from super- to 
subcritical) for a class of systems, which includes the model in |6l was given in lfl3ll . In the present paper, 
we show that the necessary ingredients for this scenario are naturally embedded in the structure of type II 
conductance based models of excitable cells. The latter covers many important biological models, including 
the Hodgkin-Huxley (HH) model (which we use below to illustrate the proposed mechanism) among many 
other conductance based models of neurons. We propose a simple geometric mechanism that explains the 
origins of chaotic dynamics, using a family of one-dimensional (ID) unimodal first return maps (Fig. [3]>. 
The qualitative features of the Poincare maps such as unimodality and the presence of the boundary layer, 
where the map is strongly expanding, follow from the proximity of the continuous system to a degenerate 
AHB and the slow-fast character of the global vector field. Therefore, the bifurcation scenario and the mech- 
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Figure 1: The time series of v (t) of the modified HH model (for the parameter values, see Appendix). In the 
present regime, the equilibrium of £T|), Q is near the AHB. The bifurcation is subcritical but is very close 
to the border, where the type of the bifurcation changes to a supercritical. The proximity of the AHB to the 
border of criticality combined with the slow-fast structure of (fl}, (f2]) results in the irregular MMOs shown 
in this figure. 

anism of chaotic oscillations described below, are independent from the specific details of the model, but are 
determined by the bifurcation structure pertinent to a whole class of models. Furthermore, we complement 
and confirm the conclusions of [13 ] and of this paper on the presence of the chaotic attractors in the cor- 
responding families of models by computing the statistics associated with the irregular MMOs. The latter 
clearly show that the timings of spikes in the mixed mode patterns are distributed geometrically in accord 
with our theoretical estimates. 



2 The model 



A large class of models of excitable cell membranes relies on the conductance based formalism due to 
Hodgkin and Huxley iflOll . Despite a rich diversity featured by different mathematical models used to elu- 
cidate the electrophysiology of various neural cells, many of them share the same structure reminiscent of 
the classical HH model. Specifically, a typical differential equation model of a (point) neuron consists of a 
current balance equation and the equations for gating variables, which often present disparate time scales. 
In identifying new dynamical structures generated by the models of this class, the specific choices of the pa- 
rameter values in a model at hand may not be so important (unless one is interested in a particular cell type), 
as they change from one model to another and also depend on many other factors (such as temperature, 
exposure to neuroactive substances, etc). However, it is critical to preserve the overall structure common 
to conductance-based models. This is the motivation and the justification for the choice of the model for 
the present study. The model to be introduced below should not be viewed as a model of a certain neuron, 
but rather as a representative example of a conductance based model. We consider a modified HH model 
introduced by Doi and Kumagai 0. This is a system of nonlinear ordinary differential equations for the 
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membrane potential, v, and two gating variables h and n: 

C m v = F(v,h,n,I ap ), (1) 
t s s = — — ,se{h,n}, (2) 

where F(v, h, n, I ap ) := -gNamoo(v) 3 h(v - E Na ) - g K n 4 (v - E K ) - g L (v - E L ) + I ap is the sum of ionic 
currents and applied current I ap . Constants E c and g c , c € {Na, K, L} stand for reversal potentials and 
maximal conductances of the corresponding ionic currents. Functions m^v), Soo(^) and t s (v), s S {h, n} 
are used in the descriptions of the ionic channels' kinetics and have typical for the HH model forms. For 
the analytic expressions of these functions and the values of parameters, we refer the reader to the Appendix 
to this paper (see also |3 ]). The only approximation used in the original HH model to obtain CQ),© is the 
steady state approximation for gating variable m sa m^v), which uses the separation of the timescales 
in the HH model to eliminate the equation for m. A mathematical justification of this approximation via a 
center manifold reduction can be found in |[T9l . Following 0, we view fh n as control parameters, thus, 
allowing for a range of timescales in (Q]) and ©. Note that ionic conductances generating electrical activity 
in different cells have widely varying time constants. Therefore, treating fh, n as control parameters allows 
one to investigate the qualitative dynamics for a class of conductance based models, while preserving the 
biophysical structure of the HH model. In particular, the modified HH model proved very useful for studying 
mechanisms for MMOs generated by conductance based models (3][T9j[20]]. The numerical simulations of 
the modified HH model EH revealed multiple parameter configurations producing chaotic MMOs. Below, 
we describe a general mechanism responsible for generating chaos in this and similar models. 



3 The Poincare map 

We start by describing the qualitative structure of the model. From (Q} and © one easily finds the equations 
for the fixed points: 

F(v, /loo (v), noo(v), I ap ) = 0, s = s 00 (v), se{h,n}. (3) 

In the range of parameters of interest, system of equations (Q~K© has a unique fixed point residing near the 
AHB. In the present study, we fix = 1 and choose I ap so that the linearization of ([T]),© around (v, h, n) 
has three eigenvalues: 

A < and a± (4) 

where < a « > 0. Thus, ©,© is close to the AHB. We will use f n to control the type of the 
AHB (sub- or supercritical). By a linear change of variables, ©,© can be rewritten as 

(5) 

where M and DM vanish at the origin. Near the origin, © has a ID local stable manifold Wf oc (0) and 2D 
local weakly unstable manifold, W z " c (0). 

To describe the mechanism for chaotic oscillations in ©, we refer to the plot of a typical phase trajectory 
near the origin. The trajectory shown in Fig. [2] approaches a weakly unstable equilibrium at the origin very 
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Figure 2: A part of the phase trajectory of © is plotted during the phase of the small oscillations between 
two successive spikes. Prior to leaving a neighborhood of the origin the trajectory remains close to a 2D 
slow manifold. The intersection of the slow manifold with cross-section E (used for the construction of the 
Poincare map) is indicated by the line of circles. For computational convenience, in the present study we 
assumed that E intersects the slow manifold along the curve of maxima of r\\. Other suitable choices of E 
are expected to yield qualitatively similar results. 

closely along W s (0) and then spirals out along W u (0). Just near the origin the oscillatory dynamics 
can be understood via the local analysis of ([5]) (cf. [13]). The oscillations of the intermediate amplitude 
shown in Fig. [2]display distinct canard structure. Note that different rings of the spiraling trajectory lie very 
closely to each other along the dashed curve. The dashed curve in Fig. [2] indicates the canard trajectory 
separating small and intermediate oscillations from large amplitude spikes lfl4ll . Once the amplitude of 
oscillations gets sufficiently large, at the end of the long excursion along the dashed curve, the trajectory 
either escapes from a small neighborhood of the origin and generates a spike, or makes another cycle of 
oscillations near the origin. The mathematical analysis of canard trajectories in M 3 and MMOs in general 
is a difficult and technical problem (see Q] [131 El dU |20l [231). The irregular MMOs such as shown in 
Fig. [2] have not been analyzed yet. We will address the mathematical analysis of this problem in the future 
work. In the present paper, we rely on the numerical observations outlined above, which suggest that near 
the origin the trajectories remain close to a 2D slow manifold. Therefore, the corresponding dynamics can 
be analyzed with a ID Poincare map. Below, we construct such map using the combination of the analytical 
and numerical techniques. 

To this end, we introduce a cross-section, E, transverse to the slow manifold, and look at the successive 
intersections of the trajectories of (f5]) with E (Fig. [2]). The Poincare map has especially convenient form in 
terms of 

x = a{ril + . (6) 

Specifically, we define 

P :x(ti)^ x(t i+1 ), i = 1,2,..., (7) 

where t{ are the times of successive crossings of the trajectory of £T|), §2} with E. The numerically computed 
maps for two values of f n are shown in Fig. [3] Note that the small values of x correspond to the large 
amplitude oscillations of ((5]). We are interested in the oscillations of small and intermediate amplitudes. 
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Figure 3: The Poincare maps plotted for two values of the control parameter f n : 43.02 (a) and 42.88 
(b). The values of f n are chosen so that the system is close to the border between the regions of sub- and 
supercritical AHB. The values of the other parameters are the same as in the Appendix. 

Therefore, in the domain of P we distinguish the region bounded away from some 0(a) neighborhood of 
the origin, /, which supports small amplitude oscillations. We call / the principal domain of P (see Fig. [3})). 
In /, P has a distinct layered structure. In the outer region, bounded away from the origin, the map is almost 
linear. In the complementary boundary layer, the map is unimodal with the slope becoming increasingly 
steep near the origin (Fig. |3^,b). The almost linear form of the map in the outer region follows from the 
local analysis of © near the origin (cf. [ 13 ]): 

P(x) = x- 2au(x + 7) + 0(a 2 ), u = 27T/T 1 , (8) 

where a and are given in ©. Note that remarkably simple form of the map in ([8]> is afforded by the 
choice of the independent variable x ©. The first Lyapunov coefficient 7 determines the type of the AHB. 
The bifurcation is subcritical (supercritical) if 7 > 0(7 < 0). The dynamics of P in the boundary layer 
corresponds to the canard like oscillations of intermediate amplitude in ©. The formal asymptotic expres- 
sion for P in the boundary layer can be found in lfl~3l . Providing a rigorous analytical description of P in 
the boundary layer is beyond the scope of this paper: it requires the information about the canard cycles 
of ©. For the purposes of the present study, we refer to the numerical results in Fig. |3l which provide a 
clear evidence in favor of the unimodality of P in the boundary layer. The latter combined with the additive 
dependence of P on 7 (see ([8]>) suggests the mechanism for chaotic MMOs in (f5]) in the regime 7 = 0(a), 
i.e. near the border of criticality of the AHB. 

Note that for negative 7 = O(l), ([7]) has a unique fixed point, which lies in the outer region (see ([8]> for 
the description of P in the outer region). This fixed point corresponds to a stable periodic orbit of © born 
at a supercritical AHB. For increasing values of 7, by ©, the graph of the map moves down and the fixed 
point looses stability via a period-doubling bifurcation as it enters the boundary layer (Fig. [3^). This triggers 
the PD cascade (not shown in this paper; see (13], where it is shown for a closely related model) and leads 
to the formation of the chaotic attractor in accord with the classical scenario known for unimodal maps Q. 
Now we are in a position to describe the mechanism for chaotic MMOs. Note that as long as min x6 / P(x) 
remains bounded away from some O(o) neighborhood of (see Fig. [3^), P(I) C / and the dynamics is 
trapped in /. This means that the amplitude of oscillations can not get too big (recall the definition of x in 
©). However, as 7 is approaching (7 = 0(a)), the graph of P(x) over / moves down and at some point 
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Figure 4: The expected value of the number of small amplitude oscillations, MNf n as a function of r n . 
Numerical data are fitted according to theoretical estimate © with a>2 ~ 2.34 and f * « 42.983. 

develops a small window of escape (see the small horizontal segment in the graph of P near the former point 
of minimum of P, Fig. [3})). Now the irregular dynamics in I can be interrupted if the trajectory hits the 
window of escape. This event, corresponding to a spike in the continuous system (Fig. [JJ), is followed by 
the reinjection of the trajectory back to I. This in turn initiates a new series of small oscillations. Note that 
just near the transition from the subthreshold oscillations to MMOs, the size of the window is very small 
and a typical trajectory stays in I for a long time before escape. Thus, just near the transition point, the 
typical patterns consist of large (random) number N fn of small amplitude oscillations between successive 
spikes. Moreover, assuming that just before the transition, the dynamics in I is mixing, one can show that 
Nf n is distributed geometrically and estimate the parameter of the geometric distribution (cf. ifTTIl "). In the 
following section, we support the proposed scenario with results of numerical simulations. 



4 The statistics of irregular MMOs 

In qualitative form, the mechanism generating irregular MMOs in (f5]) via escape from the principal domain 
was studied in IfTTIl in the context of chaotic bursting in a class of conductance based models. Note that the 
differential equation models considered in IfTTIl and in the present paper possess quite different bifurcation 
structures. The affinity between the mechanisms for chaotic dynamics in these classes of models shows 
up only after the reduction to suitable ID Poincare maps. While the basis for the reduction to the map is 
specific to each case, the qualitative properties of the resultant maps are remarkably similar. In particular, 
the analysis of the statistical properties of the irregular bursting patterns in IfTTIl carries over to that of chaotic 
MMOs. There are two main outcomes of the analysis in IfTTIl . Translated to the present problem, they can be 
formulated as follows. Suppose f% is a critical value of the control parameter, separating the subthreshold 
oscillations (f n > f°) from the MMOs (f n < f°). Suppose, in addition that just before the transition to 
MMOs the dynamics is mixing. This assumption is motivated by the presence of the steep segment in the 
graph of the map in the boundary layer near (see Fig. [3]), as well as by the results of direct simulations of 
([I]),©- Then for f% — f n > sufficiently small, the number of subthreshold oscillation between successive 
spikes, Nf n , can be approximated by a geometric random variable (cf. ifTTlO . Moreover, by estimating 
the parameter of the geometric distribution, one can predict the functional form of the dependence of the 
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Figure 5: The probability distribution of the number of oscillations between two successive spikes in a 
typical trajectory of CD,© for f n sa 42.95. The histogram for N Tn is obtained using the numerical data 
for 10 3 interspike intervals and is fitted with the exponential probability density function approximating a 
shifted geometric distribution: 7 + Geom(66.5) (plotted with solid line). 

expected value EiVr n on the distance from the transition point f% — f n > 0. In particular, the results of ifTTTl 
imply that near the border of criticality, 



for some positive constants a\ , a2 and for a certain value of f~ , which is less than critical value f° but is 
sufficiently close to it. We verified these prediction numerically for a set of values of f n near the transition 
point. In all cases, we found that Nf n are distributed geometrically. A representative plot of the probability 
density function is shown in Fig. [5] Moreover, the data in Fig. @]show a very good fit with the theoretical 
estimate in the second line of ©. We did not attempt numerical verification of the estimate in the top line of 
d9j), because it requires numerical integration of CD,© over extremely long intervals of time and is expected 
to hold in a rather narrow interval of f n . In addition to verifying (©, the numerical experiments presented in 
this section support the claim that the dynamics of ©,© near the border of criticality is chaotic. In partic- 
ular, the fact that the distribution of N fn is geometric provides a posteriori justification for the assumption 
that dynamics of P is mixing in /, which was used in iTTTTl to infer that N fn is approximately a geometric 
random variable and to derive ©. Therefore, the map-based description of the transformations in the dy- 
namics of ([]]) and ©, which accompany the approach of the border between regions of supercritical and 
subcritical AHB, combined with the numerical results presented in this section form a convincing evidence 
of the existence of the chaotic attractor in this regime. We emphasize the universal character of this effect: 
the principal ingredients necessary for its realization such as the proximity to a degenerate AHB and the 
reinjection mechanism are expected to be present in a large class of systems. In particular, the same mech- 
anism is responsible for generating chaotic MMOs in a model of solid fuel combustion lfl3l . The results of 
the present study suggest that it is relevant for an important class of biological models. 
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5 Discussion 



Substantial efforts have been made to identify and to understand the mechanisms generating chaotic dy- 
namics in physical and biological systems, and in particular, in conductance based models of neural cells 
(2j [8j [TTJ [16j UQ. In the context of neuronal modeling, understanding the sources of chaos is important 
for elucidating the origins and characterizing the properties of patterns of irregular electrical activity in the 
nervous system. The latter are quite common on both cellular and network levels. There is a compelling 
evidence showing that irregular firing in many cases results from the intrinsic properties of neural cells and 
their networks, rather than from random fluctuations. Although, there is no complete understanding for 
the functional role of chaos in neural systems, there are several lines of evidence indicating its potential 
importance. First, there are modeling studies suggesting that chaotic elements enhance the information pro- 
cessing properties of the networks that they form ifToll . Second, the loss of regularity in certain patterns of 
neural activity are known to accompany the onset of some pathological states, such as epileptic seizures and 
Parkinsonian states. Therefore, better understanding of the origins of chaos in neuronal models may help to 
predict the onset of pathological states and to control them. Single cell models such as CD,© provide the 
simplest biophysically meaningful framework for studying chaos in neuronal dynamics. Our study supports 
the view that chaotic dynamics is a likely attribute of a system operating on the edge of instability |[T6l . 
Note that on the border of criticality, the regime of interest in our study, type II models acquire remarkable 
pattern-forming capacity and flexibility: upon relatively small changes of parameters they exhibit a rich 
repertoire of qualitatively different patterns ranging from subthreshold oscillations (supercritical AHB) to 
a variety of periodic and aperiodic MMOs (subcritical AHB) (3l[T3l. Thus, the ability of the system to 
detect and to encode small changes in the stimulus is greatly enhanced in this regime. The present paper 
points out that in a broad class of type II models of neural cells, the proximity to the border of criticality 
of the underlying AHB, implies chaotic dynamics. Moreover, the results in [13 ] suggest that under certain 
natural constraints, for small a > 0, this is the only region in the parameter space of ([]]), ©, where robust 
chaotic oscillations are expected. For this regime, we provide a geometric explanation for the mechanism 
of chaotic dynamics, using a family of ID Poincare maps. We conducted a series of numerical experiments 
to verify our theoretical predictions and found that the former are in a good agreement with the latter. Our 
results can be used to locate chaotic regimes in type II models and to describe the statistical properties of the 
resultant oscillatory patterns. Finally, we emphasize that the results of the present study apply to a class of 
conductance based models (for which ©,© is on ly a representative example), as well as to models outside 
biology. In particular, we have obtained qualitatively similar results (including the statistical data as shown 
in Fig. @]and[5]> for a model of solid fuel combustion (HQS. The type of the AHB in many conductance-base 
models (including the original HH model) can be switched from sub- to supercritical by varying parameters 
in physiologically relevant ranges lfT8l . Therefore, the results of this paper are relevant to explaining firing 
patterns generated by this class of models. 

Acknowledgments. One of the authors (GM) would like to thank Victor Roytburd for useful conversations. 
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Appendix. Parameter values. 



In this appendix, we list the expressions of various functions, which appear on the right hand sides of (Q]) 
and ©: 



1 



»(") = a a (v)+/3 8 ( V y se ^' n ^' c ^ = oc^Bw c£ { m '"-^ 
0.1(25 -u) 



a n (v 



exp{^}-r 

0.01(10 -v) 
exp{i^}-r 



Hc(") 

(3 m (v) =4exp{^} 



o^h{v) = 0.07 exp 



10 

>-8) 
20 



A»(v) =0.125 exp{^} 

= exp { 4F} + l 



Here, u, the membrane potential, is in mV. The gating variables n and h are nondimensional. The values 
and the units of the parameter values in (Q} and Q, which were used to generate Fig. [I]|5]are given in the 
following table: 



lap 


0.6 /iA/cm 2 


Th 


1 


9Na 


120 mS/cm^ 


9K 


36 mS/cm 2 


9L 


0.3 mS/cm 2 


ENa 


115 mV 


E K 


-12 mV 


E L 


10.599 mV 



The time constant r ra is viewed as a control parameter. We used r n = 42.925 to generate Fig. Q]and Fig. |2] 
For the values of f n used for other plots, we refer the reader to the captions of the corresponding figures. 
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